module mesh_file

    use mesh
    
    implicit none
    
contains

    subroutine mesh_file_writer(file_name)
        character(*), intent(in) :: file_name
        
        integer i, j
        
        open(10, file=file_name, form="binary")
        ! Global mesh parameters
        write(10) num_dim, mesh_class
        write(10) is_domain_regular, is_mesh_uniform
        ! Grid
        write(10) num_grid
        do i = 1, num_cell
            write(10) grids(i)%x
        end do
        ! Cell
        write(10) num_cell
        do i = 1, num_cell
            write(10) cells(i)%num_vertex
            write(10) cells(i)%vertex_x
            write(10) cells(i)%vertex_index
            write(10) cells(i)%num_neighb_cell
            write(10) cells(i)%neighb_cell_index
            write(10) cells(i)%is_neighb_cell_far
            write(10) cells(i)%num_edge
            do j = 1, cells(i)%num_edge
                write(10) cells(i)%edges(j)%num_vertex
                write(10) cells(i)%edges(j)%local_vertex_index
                write(10) cells(i)%edges(j)%norm
            end do
        end do
        ! Boundary
        if(is_domain_regular .and. is_mesh_uniform) then
            write(10) num_grid_, num_cell_
            write(10) boundary_x, span, boundary_class, dx
        else
            write(10) num_boundary_cell
            do i = 1, num_boundary_cell
                write(10) boundary_cells(i)%num_vertex
                write(10) boundary_cells(i)%vertex_x
                write(10) boundary_cells(i)%num_neighb_cell
                write(10) boundary_cells(i)%neighb_cell_index
                write(10) boundary_cells(i)%num_edge
                do j = 1, boundary_cells(i)%num_edge
                    write(10) boundary_cells(i)%edges(j)%num_vertex
                    write(10) boundary_cells(i)%edges(j)%local_vertex_index
                    write(10) boundary_cells(i)%edges(j)%norm
                end do
                write(10) boundary_cells(i)%boundary_class
            end do
        end if
        close(10)
        
        return
    end subroutine mesh_file_writer
    
    subroutine mesh_file_reader(file_name)
        character(*), intent(in) :: file_name
        
        integer i, j
        
        open(10, file=file_name, form="binary")
        ! Global mesh parameters
        read(10) num_dim, mesh_class
        read(10) is_domain_regular, is_mesh_uniform
        ! Grid
        read(10) num_grid
        allocate(grids(num_grid))
        do i = 1, num_grid
            allocate(grids(i)%x(num_dim))
            read(10) grids(i)%x
        end do
        ! Cell
        read(10) num_cell
        allocate(cells(num_cell))
        do i = 1, num_cell
            read(10) cells(i)%num_vertex
            allocate(cells(i)%vertex_x(num_dim,cells(i)%num_vertex))
            read(10) cells(i)%vertex_x
            allocate(cells(i)%vertex_index(cells(i)%num_vertex))
            read(10) cells(i)%vertex_index
            read(10) cells(i)%num_neighb_cell
            allocate(cells(i)%neighb_cell_index(cells(i)%num_neighb_cell))
            read(10) cells(i)%neighb_cell_index
            allocate(cells(i)%is_neighb_cell_far( &
                num_dim,cells(i)%num_neighb_cell))
            read(10) cells(i)%is_neighb_cell_far
            read(10) cells(i)%num_edge
            allocate(cells(i)%edges(cells(i)%num_edge))
            do j = 1, cells(i)%num_edge
                read(10) cells(i)%edges(j)%num_vertex
                allocate(cells(i)%edges(j)%local_vertex_index( &
                    cells(i)%edges(j)%num_vertex))
                read(10) cells(i)%edges(j)%local_vertex_index
                allocate(cells(i)%edges(j)%norm(num_dim))
                read(10) cells(i)%edges(j)%norm
            end do
        end do
        ! Boundary
        if(is_domain_regular .and. is_mesh_uniform) then
            read(10) num_grid_, num_cell_
            allocate(boundary_x(2,num_dim))
            allocate(boundary_class(2,num_dim))
            allocate(dx(num_dim))
            read(10) boundary_x, span, boundary_class, dx
        else
            read(10) num_boundary_cell
            allocate(boundary_cells(num_boundary_cell))
            do i = 1, num_boundary_cell
                read(10) boundary_cells(i)%num_vertex
                allocate(boundary_cells(i)%vertex_x( &
                    num_dim,boundary_cells(i)%num_vertex))
                read(10) boundary_cells(i)%vertex_x
                read(10) boundary_cells(i)%num_neighb_cell
                allocate(boundary_cells(i)%neighb_cell_index( &
                    boundary_cells(i)%num_neighb_cell))
                read(10) boundary_cells(i)%neighb_cell_index
                read(10) boundary_cells(i)%num_edge
                allocate(boundary_cells(i)%edges(boundary_cells(i)%num_edge))
                do j = 1, boundary_cells(i)%num_edge
                    read(10) boundary_cells(i)%edges(j)%num_vertex
                    allocate(boundary_cells(i)%edges(j)%local_vertex_index( &
                        boundary_cells(i)%edges(j)%num_vertex))
                    read(10) boundary_cells(i)%edges(j)%local_vertex_index
                    allocate(boundary_cells(i)%edges(j)%norm(num_dim))
                    read(10) boundary_cells(i)%edges(j)%norm
                end do
                read(10) boundary_cells(i)%boundary_class
            end do
        end if
        close(10)
        
        return
    end subroutine mesh_file_reader

end module mesh_file